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Conservative properties of finite 
difference schemes for incompressible flow 

By Youhei Morinishi 1 

i 

f 

1. Motivation and objectives 

The purpose of this research is to construct accurate finite difference schemes for 
incompressible unsteady flow simulations such as LES (large-eddy simulation) or 
DNS (direct numerical simulation). 

Experience has shown that kinetic energy conservation of the convective terms 
is required for stable incompressible unsteady flow simulations. Arakawa (1966) 
i showed that a finite difference scheme that conserves the enstrophy in the absence 

of viscous dissipation is required for long-time integration in the two-dimensional 
vorticity-streamfunction formulation. The corresponding conserved variable is ki- 
netic energy in velocity-pressure formulation, and some energy conservative finite 
difference schemes have been developed for the Navier-Stokes equations in three di- 
mensions. Staggered grid systems are usually required to obtain physically correct 
pressure fields. The standard second order accurate finite difference scheme (Harlow 
i & Welch 1965) in a staggered grid system conserves kinetic energy and this scheme 

has proven useful for LES and DNS. However, the accuracy of the second order 
i finite difference scheme is low and fine meshes are required (Ghosal 1995). Spectral 

! methods (Canuto et al. 1988) offer supreme accuracy, but these methods are lim- 

ited to simple flow geometries. Existing fourth order accurate convective schemes 
l (A-Domis 1981, Kajishima 1994) for staggered grid systems do not conserve kinetic 

energy. Higher order staggered grid schemes that conserve kinetic energy have not 
been presented in the literature. 

J The conservation of kinetic energy is a consequence of the Navier-Stokes equations 

I for incompressible flow in the inviscid limit. In contrast, energy conservation in a 

discrete sense is not a consequence of momentum and mass conservation. It is 
possible to derive numerical schemes that conserve both mass and momentum but 
do not conserve kinetic energy. It is also possible to derive schemes that conserve 
kinetic energy even though mass or momentum conservation are violated. 

In this report, conservation properties of the continuity, momentum, and kinetic 
energy equations for incompressible flow are specified as analytical requirements for 
a proper set of discretized equations. Existing finite difference schemes in staggered 
grid systems are checked for satisfaction of the requirements. Proper higher order 
accurate finite difference schemes in a staggered grid system are then proposed. 
Plane channel flow is simulated using the proposed fourth order accurate finite 
difference scheme and the results compared with those of the second order accurate 
Harlow and Welch (1965) algorithm. 

1 Permanent address: Nagoya Institute of Technology, Japan 
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2. Accomplishments 

2.1 Analytical requirements 

The continuity and momentum equations describe the motion of incompressible 
flow. For convenience later in the analysis, these equations axe written symbolically 
as 

{Coni.) = 0 (1) 

dv ' 

+ (Conv.)i + ( Pres.)i + (Visc.)i = 0 (2) 

where 

(Cont.) = g, (Pres.), = g, (Vmc.)< = (3), (4), (5) 

Here, Vi is velocity component, p is pressure divided by density, and is viscous 
stress. Henceforth, p will be referred to as pressure. 

The conservation properties of Eqs. (1) and (2) will now be established. Note 
that Eq. (2) is in the following form. 


d<f> 

dt 


+ 1 Q+ + 2 Q* + 3 Q* + ••• =0 


( 6 ) 


The term k Q * is conservative if it can be written in divergence form 


k Q* = V-(*F*) = 


dxj 


( 7 ) 


To see that the divergence form is conservative, integrate Eq. (6) over the volume 
and make use of Gauss’s theorem for the flux terms k = 1, 2, • • *, all of which are 
assumed to satisfy Eq. (7) 

^ j J f v <t>dV = - 1 J('F* + 2 F+ + 3 F+ + ( 8 ) 


From Eq. (8), we notice that the time derivative of the sum of (j) in a volume V 
equals the sum of the flux k F * on the surface S of the volume. In particular, the 
sum of <j> never changes in periodic field if k Q * is conservative for ail k . 

Note that the pressure (Pres.)i and viscous terms ( Visc.)i are conservative a 
priori in the momentum equation since they appear in divergence form. The con- 
vective term is also conservative a priori if it is cast in divergence form. This is 
not always the case, however, and we shall investigate alternative formulations. To 
perform the analysis, we regard ( Conv.)i as a generic form of the convective term 
in the momentum equation. At least four types of convective forms have been used 
traditionally in analytical or numerical studies. These forms axe defined as follows. 
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, AJ \ - dv i 

( Adv ' )i = Vj — 


(Skew)=-^- + -v^- 

{bkew.), _ 2 Qx _ + 2 v, dx . 

. „ ( dvi dv, \ 1 dvjvj 

(Rot.)i = vj (^- - -fa.) + 2~fa~ 


( 10 ) 

( 11 ) 

( 12 ) 


As mentioned above, the divergence form, ( Div .), is conservative a priori. { Adv.)i , 
[Skcw.)i, and ( Rot.)i are referred to as advective , skew-symmetric, and rotational 
forms respectively. The four forms are connected with each other through following 
relations. 

(Adv.)i = ( Div.)i — Vi ■ ( Cont .) (13) 

{Skew.), = \{Div.)i + \{Adv.)i (14) 

(. Rot.)i = (Adv.)i (15) 


We notice that there are only two independent convective forms, and the two are 
equivalent if (Cont.) — 0 is satisfied. It is also apparent that the advective, skew- 
symmetric, and rotational forms are conservative as long as the continuity equation 
is satisfied. 

The transport equation of the square of a velocity component, tq 2 /2, is Vi times 
i = 1 component of Eq. (2). 


— + vx • (Conv.)i + * (Pres.) i + v\ • (Fisc.)i = 0 (16) 

at 

In the above equation, the convective term can be modified into the following forms 
corresponding to those in the momentum equation. 


V\ ■ {Div.) 1 = ' ( Cont. ) (17) 

vj • (Adv.)i = ~ ~ ^i 2 • {Cont.) (18) 

v\ • ( Skew.)i = (19) 


Note that the skew-symmetric form is conservative a priori in the velocity square 
equation. Since the rotational form is equivalent to advective form, the four con- 
vective forms are conservative if (Cont.) = 0 is satisfied. 

The terms involving pressure and viscous stress in Eq. (16) can be modified into 
following forms. 

/ d \ dP v i dv * 

v\ • (Pres.)i = 


dx\ 


-P 


dxi 


( 20 ) 
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Terms 

in Momentum Eq. 

T ran sport Equations 

Vi 

vi/2 

K 

(. Div .) | 

O 

0 

0 

( Adv .) = (Rot.) 

0 

0 

0 

(Skew.) 

0 

O 

0 

(Pres.) 

O 

X 

0 

(Vise.) 

O 

X 

X 


Table 1 . Conservative properties of convective, pressure, and viscous terms in 
the Uj, Vj/2, and K equations. 0 is conservative a priori , Q is conservative if 
( Cont .) = 0 is satisfied, and x is not conservative. 


v\ * ( Visc.)\ 


dr\jV\ dvi 

dxj TlJ dxj 


( 21 ) 


These terms are not conservative since they involve the pressure-strain and the 
viscous dissipation. 

We can determine the conservative properties of V 2 2 /2 and v$ 2 / 2 in the same 
manner as for v\ 2 /2. 

The transport equation of kinetic energy, K = ViV{/ 2 , is times i component of 
Eq. ( 2 ) with summation over i. 


dK 

- 5 — + ■ ( Conv.) t + ♦ (Pres.)i + Vi ■ ( Visc.)i = 0 

at 


( 22 ) 


In Eq. ( 22 ), the conservation property of the convective term is determined in the 
same manner as for V \ 2 / 2. In addition, the terms involving pressure and viscous 
stress in Eq. ( 22 ) can be modified into following forms. 


Vi • (Pres.)i 


Vi • (Visc.)i 


dpvi 
dxi p 

• {Cont.) 

(23) 

dTijVi 

dxi 

dv, 

(24) 


The pressure term in Eq. ( 22 ) is conservative if (Cont,) = 0 is satisfied. The viscous 
stress term in Eq. (22) is not conservative because the second term on the right-hand 
side of Eq. (24) is the energy dissipation. 

Table 1 provides a summary of conservative properties of convective, pressure 
and viscous terms in the transport equations of u,*, v\/2 and K for incompressible 
flow. The final goal of this work is to derive higher order accurate finite difference 
schemes that satisfy these conservative properties in a discretized sense. 


2.2 Discretized operators 

Before starting the main discussion, discretized operators need to be defined. In 
this report, the discussion of the discretized equations will be limited to uniform 
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FIGURE 1. Staggered grid system in Xi — x 2 plane. 


grid systems. The widths of the numerical grid in each direction, hi, h 2 , h 3 , are 
constant. The grid system shown in Fig. 1 will be referred to as a staggered grid 
system. In the staggered grid system, the velocity components U, ( i = 1,2,3) 
are distributed around the pressure points. The continuity equation is discretized 
centered at pressure points. The momentum equation corresponding to each velocity 
component is centered at the respective velocity point. 

Let the finite difference operator acting on 4> with respect to x\ and with stencil 
n be defined as follows. 


6 n xi 


4>{x 1 + nhi/2, x 2 , x 3 ) — <f>(x 1 -nhi/2, x 2 , x 3 ) 


*1 > * 2 , *3 


nh\ 


(25) 


Also, define an interpolation operator acting on <f> in the x x direction with stencil n 
as follows. 


-j nx i 

<P 


*1, Z2, 13 


^(xi + nhi/2, x 2 , x 3 ) + ^(xi — nhi/2, x 2 , x 3 ) 

2 


(26) 


In addition, define a special interpolation operator of the product between <j> and xj> 
in the x 3 direction with stencil n. 




nx x | 


= o^( Xl + nh i/2, * 2 , x 3 )^(xi — nhi/2, x 2 , x 3 ) 

X U *2, *3 * 

+ + nhi/2, x 2 , x 3 ) ^(xi - nhi/2, x 2 , x 3 ) 


(27) 


Equations (25) and (26) axe second order accurate approximations to first deriva- 
tive and interpolation, respectively. Combinations of the discretized operators can 
be used to make higher order accurate approximations to the first derivative and 
interpolation. For example, fourth order accurate approximations lure as follows. 


9 Si4> _ 1 6 3 <t> _ d<f> Z h 4 
8^jxj 8^3X1 dx\ 640 dxj 5 


( 28 ) 
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-i lzi ~-f Xl -<{> 

8 V 8 V v 


3 dV 


T /M 4 + 


(29) 


128 dx 1 ■ 

Discretized operators in the x<i and X3 directions are defined in the same way as for 
the X! direction. 

We define two types of conservative forms in the discretized systems. k Q * in 
Eq. (6) is ( locally ) conservative if the term can be written as 


k Q 4> = 


S 2 ( k F. 


Sixj 


+ 


S 2 xj 


,20 ) 6 3 ( k F^) 

- + - 1 + •••• 


6 3 xj 


This definition corresponds to the analytical conservative form of Eq. (7). 
globally conservative if the following relation holds in a periodic field. 

EEE ‘0*^ = 0 

X\ X 2 X 3 


(30) 
k Q* is 

(31) 


The sum that appears in Eq. (31) is taken over the period of respective direction. 
AV (= h\h 2 h 3 ) is a constant in a uniform grid system. The definition of global 
conservation corresponds to the conservation property of Eq. (8) in a periodic field. 
The condition for (local) conservation satisfies the condition for global conservation. 

2.8 Continuity and pressure term in a staggered grid system 
Now we are ready to consider our main problem. First of all, let’s examine the 
conservative property of the pressure term. As we have observed, the pressure term 
should be conservative in the transport equations of momentum and kinetic energy. 

In the staggered grid system, define the discretized continuity and pressure term 
as follows. 

( Cont . — 52) = * _* = 0 (32) 


^>1 X{ 


{Pres . _ 52), s 


(33) 


The — S2 denotes that the above approximations are second order accurate in space. 
Fourth order approximations for the continuity and pressure term in the staggered 
grid system are 

^ , c^-9SiUi 16 3 U, 

(Cont. -S 4) = - t = 0, (34) 


8 

(Pres. - 54)i = ^ lP 


8 S 3 Xi 

1 6 3 p 


, • (35) 

8611, 8 6 3 Xi 

Local kinetic energy can not be defined uniquely in staggered grid systems since the 
velocity components are defined on staggered grid points. Some sort of interpolation 
must be used in order to obtain the three components of the kinetic energy at the 
same point. The required interpolations for the pressure terms in the v\ 2 and K 
equations are 

-1Ii 6 1 U i p 1Xi 


Vi 


6ip_ 

6\X{ 


6\Xi 


— p ■ (Cont — S2), 


(36) 
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FD Schemes 
far Momentum Eq. 

Transport Equations 

Ui 

Uf/2 

K 

(Pres. — 52) 

0 

X 

Oi 

{Pres. — 54) 

0 

X 

O2 


Table 2. Conservative properties of finite difference schemes for the pressure term 
in a staggered grid system. 0 is conservative a priori , Qj is globally conservative 
if {Cont. — 52) = 0 is satisfied, O 2 is globally conservative if {Cant. — 54) = 0 is 
satisfied, and x is not conservative. 


-u llR. 

8 'fai 


lx, 


- lu, s » 


-3 x , 


9 S\ Uip 




1 hUif 




— p ■ (Cont — 54). (37) 


8 6 $Xi 8 SiXi 8 63 X i 

The following relations can be used to show global conservation unambiguously. 

££!>!£ ’=£££^ (^.-52). (38) 


X\ X2 X 3 


Xi X2 X3 


£££( 

*1 * 2 X 3 \ 



61 P 

^ix, 


1*, 



hp 

63 



= £££ [/, • (Pres. - 54), 

Xi x 2 x 3 


(39) 


Therefore, Eqs. (33) and (35) are globally conservative if the corresponding dis- 
cretized continuity equations axe satisfied. 

Table 2 shows the summary of the conservative property of the discretized pres- 
sure terms in a staggered grid system. 


2.4 Convective schemes in a staggered grid system 
As we have already mentioned, local kinetic energy K (= Z7,f7;/2) can not be 
defined uniquely in a staggered grid system. Let us assume that a term is (locally) 
conservative in the transport equation of K if the term is (locally) conservative in the 
transport equations of U \ 2 / 2, U 2 2 / 2 and U$ 2 / 2 . Since the conservative properties of 
U 2 2 /2 and U 3 2 /2 are estimated in the same manner as for U\ l / 2 , only conservative 
properties of convective schemes in the momentum and U \ 2 / 2 equations need to be 
considered. 


2.4-1 Proper second order accurate convective schemes 

Define second order accurate convective schemes in a staggered grid system as 
follows. 


{Div. — 52), 


61 U^_ 

61 Xj 


(40) 


(Adv.-S2) i = U j 1Xi Y^ i 

b\X j 

(Skew. — S2)i = h Div . — S2)i + \-(Adv. — 52)< 


(41) 


(42) 
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FD Schemes 
for Momentum Eq. 

TransportEquations 

Ui 

u'i! 2 

K 

(Div. - 52) 

0 

0 

0 

(Adv. - 52) 

0 

0 

0 

(Skew. — 52) 

0 

O 

0 


Table 3. Conservative properties of proper second order accurate convective schemes 
in a staggered grid system. 0 is conservative a priori and Q is conservative if 
( Cont . — 52) = 0 is satisfied. 


(Adv. — 52); is connected with ( Div . — 52), through the following relation. 

(Adv. - S2)i = (Div - S2)i - Ui ■ (C^t~^S2) lz ‘ (43) 


(Div. — 52), is the standard divergence form in a staggered grid system (Harlow &: 
Welch 1965). (Adv. — S2)i was proposed by Kajishima (1994). (Skew. — S2)i is 
equivalent to the scheme that was proposed by Piacsek & Williams (1970). (Div. — 
S2)i is conservative a priori in the momentum equation. The product between U i 
and (Skew. — S2)i can be rewritten as 


Ui ■ (Skew. - 52)i 


Si uJ^EE^EaIJI 

6\ X j 


(44) 


Therefore, (Skew. — S 2)i is conservative a priori in the transport equation of U\ 2 /2. 
By using Eq. (43), conservative properties of the various schemes are determined. 
Table 3 shows the conservative properties of (Div. — 52);, (Adv. — 52); and (Skew.— 
52);. These schemes are seen to be conservative provided continuity is satisfied. In 
addition, the rotational form is also conservative in light of Eq. (15). 


2. 4-2 Proposal of proper higher order accurate convective schemes 

It is of interest to derive a proper fourth order accurate convective scheme for 
a staggered grid system. Existing fourth order accurate convective schemes for 
staggered grid systems ( A-Domis 1981, Kajishima 1994) do not conserve kinetic 
energy. Here, we propose the following set of fourth order accurate convective 
schemes in a staggered grid system. 


(Div. — 54); 


9 

8 6 


- lVf“) vf’> 


1 S3 

8 S3X j 


(Adv. — S4)j 


9 

8 




SjU, 

Si xj 


llj 


1 f l7T 3iCi> \ 1 

~8\S Uj 8 Uj J S 3 xj 


(45) 


( 46 ) 
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FD Schemes 
for Momentum Eq. 

TransportEquations 

u, 

Vi! 2 

K 

(Div. — 54) 

O 

0 

0 

(Adv. - 54) 

0 

0 

0 

(Skew. — 54) 

0 

o 

o 


Table 4. Conservative properties of proper fourth order accurate convective schemes 
in a staggered grid system. 0 is conservative a priori and 0 is conservative if 
(Coni. — 54) = 0 is satisfied. 


(Skew. - S4)i = \(Div. - 54), + \(Adv. - 54), (47) 

z z 

(Div. — 54), is conservative a priori in the momentum equation. The product 
between U\ and (Skew. — 54)i can be rewritten as follows. 


JJ\ ■ (Skew. — 54)i = ^ ^ 


8 xj 

l_S 3 _ 

8 6 3 xj 




(48) 


Thus, (Skew. — 54)* is conservative a priori in the transport equation of U\ 2 /2. 
The relation between ( Adv . — 54)* and (Div. — 54)* is the following. 


(Adv. - 54)* = (Div. - 54), - U t • 


Cont . — 54) lx ‘ 


Cont . — 54) 



(49) 


This equation is a proper discrete analog Eq. (13), and (Adv. — 54)*, (Div. — 54)*, 
and (5fceu\ — 54)* are equivalent if (Cont. — 54) = 0 is satisfied. Using this relation, 
the conservative properties of the present schemes are determined. Table 4 shows 
the conservative properties of the present schemes. Comparing Table 4 with Table 
1, we see that the present schemes are a proper set of convective schemes provided 
that the continuity equation is satisfied. 

Proper higher order accurate finite difference schemes in a staggered grid system 
can be constructed in the same way as for the fourth order schemes. 

2.5 Channel flow simulation 

Numerical tests of the schemes described above are performed using plane channel 
flow. The continuity and momentum equations for incompressible viscous flow are 
solved using the proper second and fourth order accurate finite difference schemes 
in a staggered grid system using the dynamic subgrid scale model (Germano et al. 
1991). The flow is drived by a streamwise pressure gradient. A semi-implicit time 
marching algorithm is used where the diffusion terms in the wall normal direction 




Youhei Morinishi 


H 


Figure 2 . LES of plane channel flow at Re=180 by proper second and fourth order 
accurate finite difference, (a) Mean streamwise velocity; (b) Velocity fluctuations. 

Symbols: : 2nd order scheme; : 4th order scheme; • : DNS, Kim, et 

al. (1987); ; U+ = 5.5 = 2.5 log y+. 
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are treated implicitly with the Crank-Nicolson scheme and a third order Runge- 
Kutta scheme (Wray 1986) is used for all other terms. The fractional step method 
(Dukowicz & Dvinsky 1992) is used in conjunction with the Van Kan (1986) type 
of pressure term and wall boundary treatment. Periodic boundary conditions are 
imposed in the streamwise and spanwise directions. 

The subgrid-scale model is the dynamic model (Germano et al. 1991) with the 
least square technique (Lilly 1992). Averaging in homogeneous directions is used. 
Filtering is performed in the spanwise and streamwise directions. 

The spatial discretization of the second order scheme is a usual one: ( Div . — 52) 
for the convective term, (Pres. — 52) for the pressure term, and ( Cord . — 52) for 
the continuity. The corresponding Poisson’s equation of pressure is solved using a 
tri-diagonal matrix algorithm in wall normal direction with fast Fourier transforms 
(FFT) in the periodic directions. The second order accurate control volume type 
discretization is used for the viscous term. 

The spacial discretization of the fourth order scheme is as follows. The convec- 
tive term, the pressure term, and the continuity are discretized by (Div. — 54), 
(Pres. — 54), and ( Coni . — 54), respectively. The corresponding Poisson’s equa- 
tion of pressure is solved using a septa-diagonal matrix algorithm in wall normal 
direction with FFT in the periodic directions. A fourth order accurate control vol- 
ume type discretization is used for the viscous term. The subgrid scale terms are 
estimated with second order finite differences. The wall boundary condition of the 
fourth order scheme is designed to conserve mass and momentum in the wall normal 
direction in a discretized sense. 

The Reynolds number based on channel half width and wall friction velocity, Re, 
is 180. The computational box is x 2 x |7r, and the mesh contains 32 x 65 x 32 
points (streamwise, wall-normal, and spanwise respectively). 

Figure 2 shows the profiles of mean streamwise velocity and velocity fluctuations 
from the proper second and fourth order schemes. Filtered DNS data (Kim et al. 
1987) are plotted as a reference in the figures. The mean streamwise velocity profile 
from the second order scheme is shifted up in the logarithmic region. This defect of 
the second order scheme is usually observed in coarse LES (Cabot 1994). Another 
defect of the second order scheme in coarse LES is the peak value of streamwise 
velocity fluctuation is too high (Cabot 1994). These defects are improved by using 
the fourth order scheme. The computational cost of the fourth order method is 
about 1.9 times that for the second order method. 

3. Future plans 

The fourth order scheme will be tested in high Reynolds number channel flow 
to see if it has a greater advantage when the velocity fluctuations have a relatively 
larger fraction of energy near the cutoff wavenumber. 
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